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Abstract. We study the drag force on uniformly moving inclusions which interact linearly with dynamical 
free field theories commonly used to study soft condensed matter systems. Drag forces are shown to 
be nonlinear functions of the inclusion velocity and depend strongly on the field dynamics. The general 
results obtained can be used to explain drag forces in Ising systems and also predict the existence of drag 
forces on proteins in membranes due to couplings to various physical parameters of the membrane such as 
composition, phase and height fluctuations. 
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1 Introduction 

Quantum field theory explains how particles can interact 
at a distance via their coupling to a quantum field [T]. 
However interaction at a distance also occurs in classical 
systems where particles or inclusions are coupled to clas- 
sical thermal fields. For instance inclusions embedded in 
lipid membranes can interact due to their coupling with 
■ the membrane height fluctuations [2] or with the local lipid 
composition 3|. As well as effective interactions between 
inclusions, coupling to a classical field can also show up 
in the transport properties of inclusions in these fields, 
notably via drag forces which can be generated and act 
upon uniformly moving inclusions. The knowledge of drag 
forces is important as they can be used to estimate ef- 
fective transport parameters, for example diffusion con- 
stants by using the Stokes-Einstein relation. Drag forces 
have been studied in a wide variety of systems, for in- 
stance on line defects moving through liquid crystals [1] 
and on dislocations in layered structures [5] and quasi- 
crystals [BJ. As well as being present for inclusions in a 
field, drag forces are also experienced by objects outside 
but interacting with the field, for instance magnetic force 
microscope tips interacting with magnetic substrates [7]. 
In a recent letter [5] we analyzed the drag on an inclusion 
which interacts with the fluctuating field, for example a 
magnetic field at a point which moves through an Ising fer- 
romagnet. Our analysis was restricted to free scalar fields 
undergoing a general class of overdamped dissipative dy- 
namics. A number of remarkable features were found for 
the drag in this class of problems (i) the average drag force 
(/) is a nonlinear function of the velocity, in general it is 
linear for small v and is characterized by a friction coef- 



ficient A = — lim„_j.o (/) jv (ii) In systems where the free 
field theory is critical (has a diverging correlation length) 
this friction coefficient can diverge and is regularized by 
the system size (iii) At large velocities the average force 
decays to zero as (/) ~ 1/v. It was also found that nu- 
merical simulations for the drag in the Ising model could 
be well fitted using results for free fields (corresponding 
to the Gaussian approximation for the field theory of the 
ferromagnet). 

In this paper we will give an extended account of the 
results and derivations of [5] . In addition we will show how 
the divergence of the friction coefficient A can be regular- 
ized by looking at the system at a finite time after the 
inclusion starts to move, rather than in its steady state, 
and show that it diverges as a power law in time. The 
fluctuations of the force about its mean value are also 
analyzed and we show that the zero velocity fluctuations 
of the force are related to the linear friction coefficient 
via a fluctuation dissipation type relation. We also pay 
particular attention to computations of drag coefficients 
in two dimensional systems. The reason for this is that 
there has been much recent interest in the diffusion con- 
stant for proteins in lipid membranes. The first theoreti- 
cal computation of the diffusion constant of a protein in a 
lipid membrane treats the protein as a solid cylinder in a 
two dimensional incompressible fluid layer hydrodynami- 
cally coupled to the surrounding bulk fluid [5]. The drag 
force on the fluid can be computed and one finds that, via 
the Stokes Einstein relation, the diffusion constant has 
a weak logarithmic dependence on the cylinder radius a. 
Here we explore the possibility that drag may be gener- 
ated by coupling to one of the several possible physical 
fields associated with the membrane, for instance height 
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fluctuations, thickness fluctuations, composition fluctua- 
tions, local phase fluctuations. As mentioned above, in [8 
it was noted that drag forces in Ising models, which are 
clearly interacting models, could be well fitted by compu- 
tations based on free field theories. In general we have no 
explanation for this, but in this paper we have analyzed 
the drag forces in the one dimensional Ising model with 
Glauber dynamics with a (weak) point-like magnetic field 
moving at constant velocity. We find an exact expression 
for (/) which turns out to have the limiting form of model 
A dynamics for a free Gaussian theory in the continuum 
limit where the correlation length is large. 



This choice of correlation function obeys the fluctuation 
dissipation relation which ensures that the equilibrium 
measure for the field is the Gibbs-Boltzmann one. 



2.2 Specific examples and applications of the free field 
model 

Before carrying out the general calculation we will give 
some examples of the sorts of field theories, interactions 
and dynamics that one can analyze with the formalism 
that follows. We start with various choices of the operator 
A. 



2 The free field model 
2.1 Model definition 

In this section we will analyze the drag force exerted on 
an inclusion moving at constant velocity v which is lin- 
early coupled to a Gaussian or free field. We denote by 
</>(r) a scalar field on a d dimensional space. We write the 
coordinates of the system as r = (x, z) where the motion 
of the inclusion is in the z direction. The Hamiltonian for 
the system is given by 



H=\jdv 0(r)A/>(r) 



hK(j>{Q{t)) 



(1) 



where Q(t) = (0, vt) is the position of the inclusion at 
time t, A is a positive self-adjoint operator and K a linear 
operator. The instantaneous force on the inclusion in the 
direction z is given by 



d 

f = h g^ K( t>\r=Q(ty, 



(2) 



i.e. it is simply the partial derivative of the total energy 
with respect to the movement of the inclusion in the direc- 
tion z, with the field values held constant. The energetic 
formulation of the way in which the inclusion interacts 
with the field thus has the clear advantage, with respect 
to say it imposing boundary conditions, of giving an un- 
ambiguous way of computing the instantaneous force. This 
energetic approach was recently employed to compute the 
thermal Casimir force in a variety of field theories with 
dissipative dynamics of the type employed here [TU] ■ 

Note that in the above we have used the operator no- 
tation 



Av(r) 



J dr'A{v 



r>(r') 



(3) 



We will consider a general over-damped dissipative dy- 
namics for the field 4> which can be written in the general 
form 



dtf>(r) „ SH 



dt 



5(j>(r) 



(4) 



where R is a positive self-adjoint dynamical operator and 
the noise is Gaussian, white in time, with correlation func- 
tion 

( V (r,t)ri(r',t')) = 2T5(t-t')R(r-r'). (5) 



A(r) = (-V 2 +m 2 )S(r), 
A(r) = (kV 4 - aV 2 )5(r). 



(0) 
(7) 



Eq. corresponds to the Gaussian approximation for 
the Hamiltonian of a ferromagnet in the Landau theory, 
as such the field (j> can be the local magnetization, the 
local composition of a binary fluid or another local order 
parameter. The form of Eq. Q comes from the Helfrich 
Hamiltonian for a lipid bilayer |llj . where <f> represents 
the height fluctuations of the membrane about its average 
height. The term k is the bending rigidity and the term a 
is the surface tension. Still at the static level, there are a 
number of choices of the coupling of the inclusion to the 
field 6 



K(r) = S(r), 
K(r) = d • V<5(r), 
K(r) = V 2 (5(r). 



(8) 
(9) 
(10) 



The coupling ([8]) is just a localized magnetic field, that in 
Eq. ([H]) is a dipole (two fields of opposite sign close to each 
other) and Eq. (JTUJ) arises in lipid membrane physics and 
represents a coupling to the membrane curvature, tend- 
ing to induce the membrane to curve upwards or inwards. 
This sort of coupling arises for proteins whose structures 
are different in the upper and lower leafs of the mem- 
brane bilayer and thus cause the membrane to become 
locally curved. For the dynamics there are a number of 
basic models that one can consider IT2L 



R(r) 
R(r) 

R(r) 



-V 2 £(r), 
1 



(2tt)3 



r/k 



exp(ik • 
477 1 k| 



(11) 
(12) 

(13) 



The dynamical operator of Eq. (|11[) corresponds to the 
simplest form of dissipative dynamics one can write down 
for a system whose total order parameter <f> = y J v dr <p(r) 
is not conserved (also referred to as model A dynamics). 
This could apply to cases such as spin systems where <f> 
represents the local magnetization, or the local phase or- 
dering parameter in lipid systems (solid, gel, liquid). The 
operator in Eq. f|12[) corresponds to the simplest dynamics 
where the total order parameter is conserved, this would 
need to be the case for systems where <j) represented the 
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local chemical composition and where the total number 
of each type of particle is conserved (model B dynamics). 
The operator of Eq. (fT5|) is less general, but has a more 
microscopic derivation, and applies to the height fluctua- 
tions of lipid membranes driven by a surrounding fluid of 
viscosity 77, in Eq. (| 13|) it is denned via its (two dimen- 
sional) Fourier transform. 



2.3 General analysis of drag forces 

With the explicit choice of the Hamiltonian in Eq. ([TJ the 
equation of motion of the field is 



dcj){v) 
dt 



-RA(j)(r) + hRK^(r - Q(t)) + rj{r,t) (14) 



where K'(r) = K(—r). 

We now decompose the field into its average part and 
its fluctuating part (f> = (4>) + ip, these two components 
obey the evolution equations 



dt 

dip(r) 
dt 



= -RA{4>(t)) + hRK^ (r - Q(t)) (15) 
= -RAip(r) +7?(r,t). (16) 



We thus see that the mean value of the field <j> depends 
on the position of the inclusion but the fluctuations about 
this mean value are independent of the inclusion. We now 
write the inclusion position as Q(i) = (0, vt) and we write 
the mean value of the field d> as 



W>(r,i)} = 9(x,z- vt,t). 



(17) 



In the coordinate system r = (x, z' — z — vt) the equation 
for g is 



The Fourier transform of g defined as 



obeys fT3] 



5(k) = J drg(r) exp(-ik • r) 



^ - ik z vg = -RAg + hRK f . 
dt 



(18) 



(19) 



(20) 



In the steady state regime we can ignore the temporal 
derivative above and find 



.9(k) = - 



hR(k)K(-k) 



R{k)A(k) - ik z v 



(21) 



In the coordinate system moving with the inclusion the 
force is given by 

(/) = h^jKg\ r=0 = J dki*,A-(k)ff(k), (22) 



and putting Eqs. (|2~Tj) and (|2"2"]l together then yields 

W-VfaS^l*™*™. (23) 



(27r) d J R(k)A(k) - ik z v 
For small v this gives 

(/) = -A« (24) 
where the coefficient of friction A is given by 



(27r) d 



R(k)A(k) 2 



(25) 



In the case where the system is isotropic, that is K, R and 
A are functions of k — Ikl we find 



A = 



dk-. 



k 2 K{kf 



(2ir) d dJ R(k)A(k) 



(26) 



We can also analyze the case where the insertion is 
inserted at a time t = and see how the force evolves in 
time. This case is especially interesting when the corre- 
sponding steady state quantities turn out to be divergent. 
Here it is convenient to work with the Laplace transform of 
the average force defined as /(s) = J Q dt exp(—st)f(t). 
Using the fact that /(0) = we can solve Eq. fTS]) by 
Laplace transforming to give 



</(«)> 



h 2 



s(27r) d 



d ^k z R(k)K(k)K(-k) ^ 



R(k)A(k) 



ik 7 v 



The static result Eq. (|2"3"]l . when it is finite, is recovered 
from the pole at s = in Eq. (|27|) . In the limit of small 
v we can define a time dependent friction coefficient X(t) 
via (f(t)) — —X(t)v. The Laplace transform of X(t) is then 
given by 



ago 



k 2 z R(k)K{k)K(-k) 



s(2n) d J (R(k)A(k) + s) 2 
and when the system is isotropic this can be written as 



(28) 



X(s) 



k 2 R(k)K 2 (k) 
s{2ir) d dj (R(k)A(k) + s) 2 



(29) 



Finally it is interesting to ask under what conditions a 
force can be generated in a direction perpendicular to the 
direction of the insertion's uniform motion in two or more 
dimensions. The calculations above can be easily extended 
to show that the force in the direction x say is given by 



</x) = 



rfk 



ik x R(k)K(k)K(-k) 
R(k)A(k) - ik z v 



(30) 



Note that this Hall like effect can be analyzed for the forces 
on vortices in superconductors, however the evaluation of 
the force in this magnetic context requires a subtle analy- 
sis of the time dependent Ginzburg-Landau equations fTJ] . 
In our problem the interaction between the inclusions and 
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the field is via a pure potential so the evaluation of the 
corresponding forces is much more straight forward. In an 
isotropic system it is clear that = 0. The perpendic- 
ular friction coefficent is given via (f±) ~ —\±v for small 
v as 



A, = 



h 2 



(27r) d 



dk 



k x k z K(k)K(-k) 
R(k)A 2 (k) 



(31) 



An interesting example of where this perpendicular fric- 
tion coefficient can be non zero is where the interaction 
term takes a dipolar form K(k) — —id ■ k, while the other 
operators remain isotropic, in this case we find 



A± = 



2h 2 d x d z 
(2n) d 



dk^ 



k 2 k 2 

Aj „, I\j „ 



R{k)A 2 (k) 



(32) 



and thus see that it can be non zero when the dipole has 
non zero components in the direction of the motion and 
perpendicular to the motion. For a fixed dipole modu- 
lus, the magnitude of the perpendicular force is maximal 
when the dipole is orientated at 45° to the direction of 
the movement. We will demonstrate the existence of this 
rather odd perpendicular force later on in simulations of 
the two dimensional Ising model. 



2.4 Regularization of divergences in the model 



The integrals appearing in Eqs. (|23l) and (|25p may diverge. 
To be more specific, we will focus on the friction coefficient 
for an isotropic system. The divergences depend on the 
dimension d of the system and the operators A, R and K. 
For small k we will take them of the form 



A(k) ~ k 5 
R(k) ~ k p 
k(k) ~ k a 



(33) 
(34) 
(35) 



In this notation, when the field theory has a finite corre- 
lation length £ = 1/m wc thus have 6 — 0. We find that 
the integral in Eq. (|25[) is infra red divergent when d < d c 
with d c given by 



25 + p - 2a - 2. 



(36) 



We note that d c increases (i) as 5 increases, i.e. long dis- 
tance excitations cost less energy, (ii) p decreases, i.e. long 
distance modes relax more quickly (iii) when a decreases, 
i.e. when the coupling of the inclusion to the field is long 
range. In the case where the drag coefficient is infra red 
divergent it is regularized by cutting off the k integration 
at an infra red cut off k m i n = n/L where L is the linear 
system size. For d < d c we find 



A ~ L dc 



-d 



(37) 



As should be expected the divergence of the friction 
coefficient also shows up in a non-analytic behavior of the 
average drag force at small v and one can show that 



when d < d c , under the conditions p + 5 > 1 and (d c — 
d)/(p + 6 - 1) < 2. 

Finally there is another way to regularize the infrared 
divergence; we can measure the friction coefficient at a 
finite time. We expect that the friction coefficient will grow 
with the time as 

\(t) ~ 1* (39) 

and one can compute the exponent 4> using the Laplace 
transform (|29| . Making the change of variable A: = s 1 ^ p+5 - > g 
and noting that the Laplace transform of t 
to s _ ( 1+ ^) , we obtain 



' is proportional 



p + S 



(40) 



where again we assume that p + 5 > 1 . 

In general the expressions given above for the drag 
force can also exhibit ultra violet divergence which must 
be regularized. There are two possible physical length scales 
which regularize the corresponding integrals 

(i) the field theory has a natural cut off k — ir/aQ where 
ao is a length scale below which the field does not fluctuate 
or beyond which its fluctuations are strongly suppressed. 
This cut off scale can be imposed by hand and taken to 
correspond to a molecular scale, for example the lipid size 
in lipid membrane bilayers, or because the Hamiltonian 
function A has corrections at higher order in k than its 
low k form given in Eq. (1351) . 

(ii) the size of the inclusion a gives a cut off k = ir/a, 
for instance instead of having a point like magnetic field 
inclusion where K — h, one can have a Gaussian dis- 
tributed field smeared over a region of size a with K = 
hexp(—2-g-). This means that the k integration is ef- 
fectively cut off at k = w/a. For the purposes of this 
paper therefore we will take the cut off to be fc max = 
min{7r/ao, n/a}. However in most cases of interest is is 
usually a which is the larger of these two ultra violet length 
scales. 

The conclusion of this analysis is that when d < d c 
the results we obtain are dominated by the long distance 
properties of the theory and we see a diverging friction 
coefficient as £ — > 0. However if d > d c the friction co- 
efficient becomes strongly dependent on the ultra violet 
cut off, for instance on the size of the inclusion. This ultra 
violet dominated regime thus lacks the universality of the 
infra red dominated regime and we must be careful in our 
choice of model and regularization to obtain physically 
meaningful results. 



2.5 Force fluctuations 

Here we will consider the statistical properties of the fluc- 
tuations of the force about its mean value. Depending 
on the system these fluctuations may be measurable and 
could provide a method for determining some of the ef- 
fective parameters describing the system. We define the 
fluctuating component of the force as 



(38) 



Af = /- (/>■ 



(41) 
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This fluctuating component can be written in terms of the 
fluctuating component of the field ip defined in Eq. (|16[) 
and is given by 



4/ = ^U Q(t) . 



(42) 



In the steady state regime the correlation function of the 
field ip is given by 

(^(r,t)V(r',0)=C(r-r',t,t') 

= T J duzl(r-u) exp(-|t-i'|AR)(u-r'). 

(43) 

Using this we find that the correlation function for the 
force fluctuation is given by 

\ JVJ J \ )l {2jr)d J z A{k) (44) 



x exp l-\t - t'\A{k)R(k) + ik z v{t - H) 
The equal time correlation function is thus given by 

W(0 )4 » = S/^»B. ,45) 



(2tt) c 



A(k) 



and is independent of the velocity v. We also see that there 
is a fluctuation dissipation like relation relating the zero 
velocity force fluctuations to the linear friction coefficient: 



dt(Af(t)Af(0))\ v=Q =T\ 



(46) 



In general any measurement of a force will not be instan- 
taneous and will depend on the temporal resolution of the 
experimental set up and will thus represent force aver- 
aged over a characteristic time scale T m associated with 
the force measurement apparatus. We define the tempo- 
rally averaged force over the time window T m as 



1 



1 7 



(47) 



clearly we have (f m ) — (/) and the variance at zero ve- 
locity G 2 m = ((f m - (f)) 2 )\ v=0 is given, for large T m , by 



2TA 
~~T ' 

± rn 



(48) 



3 Numerical simulations of drag in the Ising 
model 

In this section we perform numerical simulations of drag 
forces on inclusions in the Ising model. This is an example 
of an interacting theory where the drag forces predicted 
in free field theories should also occur. Our simulations 
in fact show that, despite their approximative nature in 
the context of interacting theories, our results for the free 
Gaussian ferromagnet account well for the phenomenology 



of drag observed in Ising systems. We will consider the 
model on a d-dimensional cubic lattice of spacing ao, with 
periodic boundary conditions, and denote by N the total 
number of sites and spins. The Hamiltonian is given by 



H = -J^S i S j -hY J K i 



> Si 



(49) 



(id) 



where J > is a ferromagnetic coupling between nearest 
neighbour spins and where hKi_i is the local field at site 
i due to the inclusion whose position is denoted by Iq. 
Here the vector Ki is the discrete version of the operator 
K of Section 2. In what follows, in order to fully inves- 
tigate the various models discussed in the paper, but to 
keep to a reasonable length and minimize the amount of 
computation time, we will restrict our study to one and 
two dimensions. 

The system dynamics is defined in the following man- 
ner: N elementary evolutions are performed during one 
unit of time. An elementary evolution consists of: 

— choosing a spin set (the way of choosing it depends on 
the dynamics and will be given below), 

— computing the energy change AH associated with flip- 
ping this set of spins, 

— flipping the spins with probability 

Pf = 1/(1 + exp(AH/T)) or leaving them unchanged 
with probability 1 — pf . 

We simulate two types of dynamics: 

— Non-conserved dynamics: only one spin is chosen ran- 
domly at each step; thus the total magnetization is 
not conserved. This choice is referred to as Glauber 
dynamics. 

— Conserved dynamics: at each step, two spins of oppo- 
site sign are randomly chosen; the total magnetization 
is conserved. This is a form of Kawasaki dynamics. 

The inclusion moves in the z direction with velocity v, 
so io,z(t) = mt(vt/ao) (int denoting the integer part) i.e. 
it performs one step every ao/v units of time. To measure 
the force in the z direction in a given configuration of the 
spins system, we compute the energy H + if the inclusion 
was at the position io + z and H- if it was at the position 
«o — z, where z is the lattice link vector in the direction 
z. The instantaneous force is then / = — H + )/2ao. 
Explicitly 



2a 



-(i +z) 



K, 



(io- 



(50) 



As we are interested in the average force, we average over 
all Monte-Carlo (MC) time steps after first achieving a 
steady state in the simulation: 



Jmc-*' 



1 



T. 



MC 



Tmc 



/(*) 



(51) 



We will present two kinds of plots of the results of 
our simulations (i) the average magnetization in the rest 
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frame of the inclusion to see how the inclusion polarizes 
the spins around it and how this polarization cloud is de- 
formed by the inclusion's movement (ii) the average drag 
force as a function of the velocity (i.e. (f)(v)). Note that 
the polarization induced by the inclusion is basically a 
generalization of the polaron in solid state physics; the 
polaron is the response of a body's polarization field due 
to the presence of an electron and modifies the dynamical 
properties of the electron [TS] . 

3.1 Point like magnetic fields in one and two 
dimensions 

Here we study the case where the inclusion creates a point 
like magnetic field: -fQ-i D = 5i^ . We take h < 0, so that 
the average magnetization is (positively) proportional to 
the potential seen by the inclusion. 

Figure [T] shows the magnetization profile for Glauber 
and Kawasaki dynamics, for four different speeds (with 
parameters (3 — 1, J — 1 and h = —10). When the parti- 
cle is at rest, it sees a spherically symmetric potential, thus 
is does not experience any net force. As the velocity in- 
creases, the profile becomes asymmetric and its amplitude 
decreases: the system has less time to react to the pres- 
ence of the inclusion and thus the polaron is deformed and 
becomes weaker. The main differences between the two dy- 
namics are (i) the magnetization far from the particle is 
not zero with Kawasaki dynamics, because the total mag- 
netization must remain zero; (ii) the polaron deformation 
appears to be larger with Kawasaki dynamics. The mean 
force is plotted against the velocity in Fig. [21 this figure 
shows that the force has a linear dependence on v for small 
v, reaches a maximum and then decreases as 1/v for large 
v. This behavior is in agreement with our general results 
for free fields: the asymmetric profile is responsible for 
the force and for small v the asymmetry increases with v 
whereas for large v the profile amplitude diminishes. As 
the deformation is larger for Kawasaki dynamics, the force 
is larger. The fits with our analytical results for model A 
and B dynamics are performed by varying three parame- 
ters: the cut-off and a dilatation for each axis. 

The 2d simulations give similar results (Figs. [21 SI). 
In two dimensions we are in the high temperature regime 
before the ferromagnetic phase transition (J3 = 1, J = 0.4, 
h = -6.66). 

3.2 Dipoles in two dimensions 

Here, the inclusion interacts as a dipole: JQ_i = (<5i.; — 
<5i,i -u), where u is a unit vector giving the direction of the 
dipole. Fig. [5] shows the drag force for dipoles perpendicu- 
lar (u = x) and parallel (u = z) to the direction of motion 
for Glauber dynamics. We also show the contour plots for 
the local magnetization profile for both cases, at velocities 
v = and v = 0.5. As seen from the force curves, at slow 
speed the force does not depend on the orientation, when 
the speed increases the force for the perpendicular dipole 
becomes larger. 



Finally, we compare the forces parallel and perpendic- 
ular to the motion for a dipole orientated at 45° to the 
direction of the motion (u = x + z) in Fig. [51 The force 
in the direction x is calculated in the same way as that 
described above for the force in the direction z. We see 
that the transverse force has the same order of magnitude 
as the longitudinal force, and the same general form. Also 
shown on the right is the corresponding contour plot of the 
local magnetization generated by the dipole at rest and 
for v — 0.5. Whereas at v = the magnetization profile 
appears antisymmetric about the direction of the dipole, 
when the dipole moves it experiences a force which pushes 
it to the left on the bottom right figure of Fig. [5] This is be- 
cause the leading component of the dipole barely sees the 
polarization created by the lower component, whereas the 
magnetization created by the leading component pushes 
away the lower component. 



4 The one dimensional Ising model with 
Glauber dynamics 

In the Section 2. we have analyzed the drag on inclusions 
in free fields. Our simulations in Section 3. were however 
for the Ising ferromagnet. We showed that the force mea- 
sured in these simulations could be remarkably well fitted 
by a free field theory. Here we show that the drag for 
a point like magnetic field in the one dimensional Ising 
model with Glauber dynamics [16] can be exactly solved 
within the linear response regime, where /3h <C 1, and that 
the force so obtained is exactly of the form predicted from 
the model A dynamics for the Gaussian ferromagnet. 

As in our simulations we will compute the symmetrized 
instantaneous force given by 

(/(«)) - ~ [(S i0 (t)+i) - (S ioi t)-i)} , (52) 

where io(t) = mt(vt) and we have set the lattice spacing 
do = 1. The time dependent magnetic field in this problem 
can be written as 

hk(t) = hoS ktio (t). (53) 

We will work in the regime where the applied field is small 
and apply linear response theory to write 

(s j (t)) = (SMo + h^f_js (gg) o W 

(54) 

where (-)o indicates averaging in the absence of the field 
h. We also assume that the dynamics of the system in 
absence of the field h evolves from a statistically homoge- 
neous initial state such that 

(5i(t)) - (Sj(t)) , (55) 
for all i and j. Along with Eq. ([541 in Eq. (1521 this yields 
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Fig. 2. Dashed lines: average drag force / in the Id Ising model as a function of v for Glauber (a) and Kawasaki dynamics (b). 
Solid lines are the fits of model A (a) and model B (b) dynamics for the Gaussian ferromagnet. 
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v = 0.25 





15 25 35 



Fig. 3. Contour plot (color online) of the magnetization profile (polaron) for the 2d Ising model with Glauber dynamics about 
a local magnetic field at a single point moving with velocity v. (high temperature phase: /3 = 1, J = 0.4, h — —6.66). 



(a) 



(b) 





Fig. 4. Dashed lines: average drag force / in the 2d Ising model as a function of v for Glauber (a) and Kawasaki dynamics (b). 
Solid lines are the fits of model A (a) and model B (b) dynamics for the Gaussian ferromagnet (/3 = 1, J = 0.4, h — —6.66). 



(/(<)> 

- H. v 

1 k J - 



ds 



SS„ 



»o(*) + l \ _ / u ^o(t)-l 

Shk(s) I \ Sh k (s) , (L 

^i (*) + l \ / 



SS„ 



8h 



*o(«) 



'«o(t)-lJ 



8h 



io(s) 



J k,i (t) 



(56) 



In addition, for a system in equilibrium, one has the fluc- 
tuation dissipation theorem |17j 

T,/. ■ , \ s dC(i,j,t, s) dC(i-j,t-s) 



(9.s 

for t > s, and where 

C(i,j,t,s) = {S i (t)S j (8)) 



ds 



(60) 



(61) 



The response function for the unperturbed system is is the spin-spin correlation function. Therefore thermal 



defined by 



equilibrium we find 



(57) 



</(*)> = 



Ph 2 



and for a system in thermal equilibrium we may write 

Tl(iJ > t,s)=n(i-j > t-s) ) (58) 

as we have spatial and time translation invariance and 
thus 



d 

ds— (C(i (t)-i (s) + l,r) 

OT 



(/(«)> 



h 2 ft 



fl^ds [1Z(i {t)-io(s) + l,t-s) 

-K(i (t)-i (s)-l,t-s)} 



(59) 



-C(io(t) - i (s) - l,r)) T=t _ s . 

(62) 

The correlation function obeys the equation 
-C(i, t) = -C(i, t) + ^ (C(i + 1, t) + C(i - 1, t)) • (63) 

where 7 = tanh(2/3J) [IB] . We now consider a system with 
2L + 1 spins at sites — L, . . . , • • • , L and periodic bound- 
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Fig. 6. Forces for a dipole orientated at 45° to the motion with Glauber dynamics. Solid line: average force perpendicular to 
the motion; dashed line: force parallel to the motion (/? = 1, J = 0.4, h = —6.66). Contour plot (color online) for two values of 
v. 
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ary conditions. We define the discrete Fourier transform 
of C via 



C(i)= E CWexp 



/ 2irijk 
\2L + l) ' 



(64) 



fc=-L 

and thus the Fourier coefficients are given by 

L 



C(k) 



1 



2L 



exp 



2iTijk 
2L+1 



(65) 



The initial condition for Eq. (|63p is given by the equilib- 
rium correlation function 



=1,1*1, 



(66) 



where 77 = tanh(/3J). We can now use Eq. (l63l) to express 
Eq. in terms of its Fourier representation to find 



(/(*)> = -M 2 

1 — 7 cos 



/ dsJ2C(k,(t-s)) 
J— 00 ,, 



/ 2?rfc 



x sm 



V2L+ 1 

/ 2irk \ /2mk(i Q (i) - i (s)) 
\2LTi J 6XP V 2L + 1 



(67) 



Now in the continuum limit where £ ;g> 1 we write simply 
that io(t) = vt and use the solution 



C(k,t) = C(k,0)exp[-t l- 7 cos 
to obtain 

C(k,0) l-7cosf 



/ 2nk 
V2L + 1 



(68) 



(f(t))=(3ih 2 '£- 



/ 2nk 
2L + 1 



sm 



/ 27rfc \ 
\2L+1 J 



l-7cos(^ 



2-nikv 
2L + 1 



(69) 

The initial condtion Eq. (|66| along with (1651) then gives 

L 



C(k,0) = 



1 



2L 



77'-" exp 



2m jk 
2L+1 



1 



1 - rj 2 



2L + 1 1 + rf - 2?ycos 



2?rfc 
2L+1 

I I 1 



2L + 1 cosh(2/3 J) 1 _ 7C0S ^^ife_ 



(70) 



where we have taken the limit of large L and assumed 
that 77 < 1 (i.e. non zero temperature). Putting this all 
together then yields 



</(*)> 



f3ih 2 



cosh(2f3J)(2L 



sm 



2irk 
2L + 1 



7 COS 



, / 2irk \ 
' ^2L+lJ 



2irikv 
2L + 1 

(71) 
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Fig. 7. Drag force in one dimension for small field in an Ising 
model with Glauber dynamics (crosses) versus analytical result 
([72|l (solid line); /3 = 1, J = 1.5, h = 0.2. 

Now the sum can be written as an integral for L large to 
give 



(/(*)> = 



(3ih 2 



2?rcosh(2^J) 



dk 



sin (k) 



1 — 7 cos (k) — iki 



(72) 



We can recover the link with the continuum models stud- 
ied here if we consider the limit where the inverse corre- 
lation length mCl. Here we have 



m = — 111(77) ~ 2 exp(— 2/3 J). 



(73) 



If we take m small the integral in Eq. (|T2")) is dominated 
by its behavior at small k, in addition we have 7^1 — 
2exp(— 4/3 J) = 1 — m 2 /2 which yields 



(/(*)> = 



j3imh 2 



dk 



k 



2ikv 



(74) 



which has the same form as that for model A dynamics 
for a Gaussian ferromagnet. 

The analytical result (|T2")) has been compared with a 
simulation; the results given Fig. [7] show a good agreement 
between them. We interpret the fact that the analytical 
result overestimates the simulations result as a trace of a 
non linear response in the field h. 



5 Application to proteins in lipid membrane 
5.1 General analysis 

Here we will try and investigate some possible sources of 
drag in lipid membrane models. A way of computing the 
diffusion constant of an insertion, such as a protein in a 
lipid membrane, is via the Stokes Einstein relation 



D 



k n T 



(75) 
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where At is the total friction on the protein. There are a 
number of possible sources of drag in lipid membranes. 
The first treatment of this problem was by Saffmann and 
Delbriick (SD) [5] who computed the hydrodynamic drag 
by treating the low Reynolds number Navier Stokes equa- 
tions for a slab of flat 2d fluid containing a solid cylindrical 
insertion. The movement of the fluid sets up a hydrody- 
namic flow and the resulting friction on the cylinder is 
computed using the stress tensor. The coupling to the 
bulk external fluid is very important and is essential to 
find a finite result for the drag, as a purely two dimen- 
sional treatment gives a divergent result coming from the 
long range nature of hydrodynamic interactions in two 
dimensions. The hydrodynamic drag computed by SD is 
given by 

Ahydro — 7 — , — ~, (.10) 



In 



7 



where a is the cylinder radius and r\ m and r\ w are the the 
viscosities of the membrane and surrounding fluids respec- 
tively. The term h represents the height of the cylinder or 
membrane and 7 ~ 0.5772 is Euler's constant. This for- 
mula is valid in the regime where a 3> h, i.e. for proteins 
which are large relative to the membrane thickness and 
when r) m ^> r) w . The coupling of the 2d flow to the 3d fluid 
is very important in this hydrodynamic treatment, for ex- 
ample if there is a hard wall in the proximity of the fluid 
membrane, the behavior of the diffusion constant changes 
to D ~ 1/a 2 [IE]. 

Recently in |19] a detailed experimental study, and 
comparison of other results in the literature, of protein dif- 
fusion constants seems to suggest that for membrane pro- 
teins and peptides the diffusion constant scales as D ~ 1/a 
(which is consequently a much stronger dependence on 
the protein radius than D ~ ln(l/a) predicted by SD). 
In [19] it is suggested that the apparent failure of the SD 
formula may be due to the fact that the membrane is 
quite heterogeneous on small length scales and that the 
model of a perfect incompressible fluid is perhaps not well 
adapted for small inclusions. It is also pointed out that 
on larger length scales thermal fluctuations and undula- 
tions may dissipate velocity gradients. Indeed extensive 
numerical simulations have shown that the coupling of the 
protein position to local membrane curvature (and hence 
height fluctuations) reduces the diffusion constant of in- 
clusions PHI and scaling arguments show that proteins 
whose hydrophobic cores are mismatched with the equi- 
librium thickness of the lipid bilayer also experience ad- 
ditional drag forces [3T] . The effect of mismatch is clearly 
seen in some of the experimental studies reported in |19] . 

In the spirit of the comments of [TH] and the study of 
|21] , we will tentatively examine various scenarios leading 
to drag forces on membrane inclusions which are linearly 
coupled to physical fields in the membrane. We should 
bear in mind the limitations of this approach. First it is 
clear that we are ignoring possible hydrodynamic flows 
created by the movement of the inclusion. If substantial 
hydrodynamics flows are established by protein movement 
then the order parameter field (depending on its precise 
nature) can be expected to be convected with the flow. 



Clearly this effect is ignored in our study, however for 
small inclusions where the inclusion pushes past the lo- 
cal lipids, rather than entraining a hydrodynamic flow, 
we expect this approximation to be good, at least con- 
cerning orders of magnitude estimates. Secondly there is 
also a nonlinear coupling between inclusion position and 
the fluctuating field (due to the absence of the field in the 
region of the inclusion) . However if the linear coupling is 
sufficiently strong the contribution of the mean field like 
term we compute here should dominate the drag due to 
nonlinear terms. 

In d = 2, the friction coefficient (|2"6]l reads, 



2?r 



(Xt\i ~ 



R(k)A 2 (k) 



(77) 



where a c is the short distance cut-off scale corresponding 
to the larger of the two scales a, the insertion size, and 
ao, the underlying cut off for the field fluctuations. In our 
models, the operators will be of the forms of Eqs. (|3"3T 
1351) : in order to be more precise, we write on dimensional 
grounds 



A(k) = n A a s 'k s ' (k 2 + to 2 ) , 
R(k)A{k) = To X a p +5 ' +2 k p+5 ' \k 2 + m 2 ) , 
K(k) = ri K a%k a , 



(78) 
(79) 
(80) 



where 11 a and iik are energies and To is a microscopic time 
scale associated with the dynamics. Here 6' is simply the 
exponent S of the previous sections when ra/0. 

Using these expressions in the integral above, and ex- 
tracting the mass dependence by setting k — mq, we ob- 
tain for the friction coefficient 

A = ^hh 2 (ma r^'g x (H^S.) , (81) 



where the function g\ is defined by 



dq 



q 



3+2a-p-2<5' 



(q 2 + I) 2 



(82) 



In the following, we assume that this integral is not in- 
frared divergent, i.e. that 26' + p — 2a < 4; this reads 
d = 2 > d c with the critical dimension (pJrJ)) . 

It now remains to determine how the amplitude of the 
interaction h should be computed. It is clear that the value 
of h should depend on the value of the size of the inclusion. 
A simple, semi-macroscopic, way of doing this, proposed 
in [5] , is the following. The energy of interaction between 
the field <j) and the inclusion is easily computed from the 
free field theory and is given by 



dk 



kk 2 {k) 



(83) 



2(2tt) Jo A(k) 

We now expect that e is a function of a and that for small 
a 

e(a) ~ ~2TTjja — iraia 2 (84) 
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where 7/ and 07 are effective (negative) line and surface 
tensions for the inclusion in the membrane due to the in- 
teraction with the field <j>. Now if we assume that a is small 
we can neglect the surface tension term and equating Eq. 
and Eq. (J83I) we find 



27T7/a 



h 2 p, 2 K (a m) 



2a-S' 



where 



2(2tt) 

9e(x) = 



dq 



1 



(86) 



Using the resulting expression for h in terms of 7/ and a 
we then obtain 



47r 7 jar ogA (^) 
al(ma Q )P+6'g e (^ 



(87) 



Before examining a number of models of proteins in- 
serted into membranes we will explore a few general con- 
sequences of the above expression. Given that we expect a 
and ao to be small we should consider the functions g\ and 
g e in the limit x — > 0. The integrals defining these func- 
tions are finite (and thus independent) of a in the cases 
where 2a — p — 25' < and 2a — 5' < respectively. In 
this case we find the generic behavior 



47r7jaroff A (0) 
A ™ a 2 (ma )P+S'g e (0)' 



In this scenario we see that the dependence of the friction 
coefficient on the inclusion size is always linear and it has 
a strong dependence on the correlation length of the field 
4>. The scaling of the friction coefficient with the inclusion 
size is A ~ a, if this drag dominates all other sources of 
drag application of the Stokes-Einstein relation gives 



D 



(89) 



Note that we will also recover this dependence on a if it is 
the case that a < a , i.e. the underlying cut-off of the field 
<f> is greater than the one corresponding to the inclusion 
size. Note that apart from the hydrodynamic case where 
p = — 1, p is positive or zero, therefore if 2a — 5' < 0, then 
2a — p — 25' < also. In the case where 2a — p — 25' > 
and which in most cases will also imply that 2a — 5' > 0, 
the integrals defining both g\(x) and g t (x) will diverge 
and we find 



9x{x) 

9e{x) 



1 



1 



2a- p- 25' x 2a -P- 2S ' 
1 1 



2a - 6' x 2a ~ s ' 



and thus 



2a — p — 25' \ irao J 



An^/iaTQ 2a — 5' 



P+S' 



(91) 



(92) 



Again if this drag dominates, it gives a diffusion coefficient 
scaling with inclusion size as 



D 



1 



aP+S'+i ■ 



(93) 



In this case we see a different dependence on the inclusion 
size a, the physics of the problem is controlled by short 
distance behavior and the drag is independent of the cor- 
relation length ( = 1/m of the fluctuating field. A final 
possible case is where 2a — p — 26' < but 2a — 5' > 0, in 
which case we find 



A ! 



47r77aT g A (0)(ma) 2a -' 5 ' 2a - 6' 
al(ma )P+ 6 ' n 2a-8' 



(94) 



this is a particularly interesting case as the friction coef- 
ficient has a strong dependence on both the correlation 
length of the field <j) and on the inclusion size. Again if 
this drag dominates the diffusion constant scaling with 
inclusion size is 

D - ^TT^ ( 95 ) 
We now discuss some specific models. 



5.2 Insertion with curvature coupled to membrane 
height fluctuations 

In [2D] the authors numerically investigated the diffusion 
constant of inclusions in model membrane systems where 
the inclusion tends to impose a preferred local curvature 
on the membrane. In their model a quadratic coupling was 
also considered. Here we consider the model with a simple 
linear coupling. The inclusion is coupled to the membrane 
height fluctuations, the Gaussian Hamiltonian of which A 
is given by Eq. ([7]), the dynamical operator R is given 
by Eq. ((13} and K by Eq. fTD]), After Fourier transform- 
ing, we obtain A(k) — nk 2 (k 2 + m 2 ) (with m = s/oJk), 
R(k) — (4:7]k)~ 1 and K(k) = —k 2 . In our general notation 
we have 5' = 2, a = 2 and p = —1 and we are in the case 
where 2a — 5' — 2 and 2a — p — 25' = 1. In terms of the 
physical parameters of this model the drag coefficient for 
small insertion size a is thus given by 



A = 



32777/a 2 



(90) an d the dominance of this drag would imply that 



D 



(96) 



(97) 



We thus see that this result is quite insensitive to the 
correlation length of the height fluctuations (and thus the 
surface tension) assuming that they are large with respect 
to the insertion size. 
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5.3 Insertion coupled to a non-conserved order 
parameter 

An insertion such as a protein or peptide can couple to 
various physical fields in a lipid membrane other than the 
height fluctuations. For instance in a lipid monolayer lo- 
cal tilt angles of the lipid heads or tails may be changed 
by the presence of an inclusion. Also if there are several 
lipid phases such as liquid gel and solid, the inclusion may 
prefer to be in one of these phases. This general idea can 
be modelled by assuming that the inclusion couples lin- 
early to the order parameter representing one of these 
fields. The simplest Hamiltonian for this order parameter 
has the form of that for the Gaussian ferromagnet where 
A(k) oc k 2 + m 2 . The simplest diffusive dynamics is given 
by model A dynamics with R{k) oc 1 and a linear coupling 
to the field 4> gives K(k) = 1. This dynamical model does 
not conserve the integrated field as there is no reason that 
it should be conserved. Here, in our general notation, we 
have 5' = 0, a = and p = and we are thus in the 
case where 2a — 5' = and 2a — p — 25' = 0. We see that 
we are in the marginal case for both functions g\ and g e . 
Furthermore we can identify the time scale tq using the 
diffusion constant for the dynamics of the field </>, Dq, via 
DqTq = a§, where ao is the lipid size and Do can be es- 
timated from the lipid translational diffusion constant or 
lipid rotational diffusion constant, depending on field in 
question (for instance if the field in question describes the 
orientational order of the lipids, then the lipid rotational 
diffusion constant could be used to give the appropriate 
time scale). Here for small x we find 



g\(x) w -ln(x) and g e (x) » — ln(x) 



which leads to 



"Do"' 



(98) 



(99) 



and gives an estimation of the insertion diffusion constant: 



D 



k B TD 



(100) 



5.4 Insertion coupled to a conserved order parameter 

The insertion may also be coupled to a conserved field, 
describing, for example, the local lipid composition in the 
case where there are several lipid types. We take the same 
A and K operator to describe the energy of the order 
parameter describing local chemical composition. However 
we now use a conserved dynamics: R is set by l|12p. giving 
R(k) oc k 2 . We thus have 5' = 0, a = and p = 2, which 
gives 2a — 5' = and 2a — p — 25' — — 2. The function g e 
is unchanged but we find g\ ~ 1/2 which gives 



A = 



Doa 2 m 2 H±) 



(101) 



where we have again used DqTq — Oq, and where -Do can be 
estimated from the lipid translational diffusion constant. 



This leads to the estimate 

k B TD Q a 2 m 2 \n(±) 



D 



2-K^fia 



(102) 



for the protein diffusion constant. We should note that 
even though there is a logarithmic correction, we would 
expect to experimentally measure D ~ \/a as the loga- 
rithmic term would require decades of length scales (thus 
leaving the realm of validity of the calculation) to detect. 
Interestingly here, in contrast with the previous models, 
we should see a strong dependence of D on the correlation 
length of the field. 



6 Conclusions 

We have seen that inclusions which are linearly coupled 
to classical fields with dissipative dynamics are subject 
to drag forces which exhibit a rather rich behavior. No- 
tably the drag force is a non monotonic function of the 
inclusion velocity v. Generically the force is a linear func- 
tion of the velocity at small v and is characterized by a 
friction coefficient A. The force then attains a maximum 
value, and for large v decays as 1/v. The force is phys- 
ically generated by the deformation of the polarization 
profile of the field about the inclusion by the inclusions 
motion. This phenomena is analogous to the way in which 
electron dynamics is renormalized by their associated po- 
laron in solid state physics |15) . The linear coefficient of 
friction A is of particular importance as it can be used to 
estimate diffusion constants via the Stokes-Einstein rela- 
tion and because it can exhibit divergent behavior when 
the corresponding field theory is critical, i.e. has a diverg- 
ing correlation length, for example at a critical demixing 
transition. 

The results we have presented are valid for free or 
Gaussian field theories. We were able to compute the drag 
for the one dimensional Ising model for a weak inclusion- 
field interaction but it would be interesting to go beyond 
the Gaussian approximation to understand the physics of 
drag forces for general interacting field theories. Having 
said this we note that the Gaussian model does seem to 
capture most of the phenomenology seen in our simula- 
tions of the one and two dimensional Ising model. 

Finally the experimental measurement of these drag 
forces presents an interesting challenge. It may be possible 
to carry out experiments using atomic force microscopy or 
magnetic force microscopy if the interaction between the 
microscopic tip and the surface can be sufficiently well 
characterized. Also, the thermal or critical Casimir force 
predicted by Fisher and de Gennes [32] has recently been 
successfully measured [2 3) in a binary fluid mixture at 
criticality. It may be that the technical advances made 
to carry out this measurement, the chemical treatment 
to tune the interaction between the fluid components and 
surfaces, and the optical force measurements could be ap- 
plied to the study of the drag problem. We note also that, 
beyond measurements of the average force, it would also 
be interesting if the predictions made here about force 
fluctuations could be verified experimentally. 
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